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Abstract 



We present a new third-order central scheme for approximating solutions of 
systems of conservation laws in one and two space dimensions. In the spirit 
of Godunov-type schemes, our method is based on reconstructing a piecewise- 
polynomial interpolant from cell-averages which is then advanced exactly in time. 

In the reconstruction step, we introduce a new third-order, compact, CWENO 
reconstruction, which is written as a convex combination of interpolants based 
on different stencils. The heart of the matter is that one of these interpolants is 
taken as an arbitrary quadratic polynomial and the weights of the convex combi- 
nation are set as to obtain third-order accuracy in smooth regions. The embedded 
mechanism in the WENO-like schemes guarantees that in regions with discontinu- 
ities or large gradients, there is an automatic switch to a one-sided second-order 
reconstruction, which prevents the creation of spurious oscillations. 

In the one-dimensional case, our new third order scheme is based on an ex- 
tremely compact four point stencil. Analogous compactness is retained in more 
space dimensions. The accuracy, robustness and high-resolution properties of our 
scheme are demonstrated in a variety of one and two dimensional problems. 

Key words. Hyperbolic systems, central difference schemes, high-order accuracy, non- 
oscillatory schemes, WENO reconstruction, CWENO reconstruction. 
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1 Introduction 

We are concerned with multidimensional systems of hyperbolic conservation laws of the 
form 
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Methods for approximating solutions to equation ( |1.1|) have attracted a lot of atten- 
tion in recent years (see ||], [1C], [24] and the references therein). 

In this work we focus on Godunov-type schemes, where one first reconstructs a 
piecewise-polynomial interpolant which is then advanced exactly in time according to 
( [Lip and finally projected on its cell-averages. Generally, one can divide Godunov-type 
schemes into two sub-classes - upwind methods and central methods. 

In upwind schemes, one first reconstructs a polynomial in every cell, which is then 
used to compute a new cell average in the same location in the next time step. This 
procedure requires solving Riemann problems at the discontinuous interfaces. For high- 
order methods, instead of analytically solving the resulting Riemann problems, one 
typically implements approximate Riemann solvers or some form of flux splitting. For 
systems of conservation laws, or in the demanding context of more space dimensions, 
this procedure turns out to be more intricate as such Riemann solvers do not exist. The 
essentially non-oscillatory (ENO) methods by Harten are the prototype of high-order 
methods (see ||, |Z| and the references therein). A recent review of ENO and WENO 



methods can be found in |21 



Central schemes, on the other hand, are based on averaging over the Riemann fans, 
a procedure which is typically done by staggering between two grids. They require 
no Riemann solvers, no projection along characteristic directions and no flux splitting. 
Therefore, all that one has to do in order to solve a problem is to supply the flux 
function. They are more simple when compared with upwind schemes. 

The major difference between different central methods is in the reconstruction step, 
where one computes a piecewise-polynomial interpolant from the previously computed 
cell- averages. 

The prototype of central schemes is the Lax-Friedrichs |1 scheme which is based on 
a piecewise-constant interpolant. Even though it is very robust, it is only first-order 
accurate and, moreover, it suffers from excessive numerical dissipation. A second-order 
central method was proposed by Nessyahu and Tadmor in |19| . This method is based 
on a MUSCL-like |J piecewise linear interpolant and nonlinear limiters which prevent 
spurious oscillations (see |2(J for a different approach). A variety of extensions to these 
methods were suggested. The one-dimensional third-order method of Liu and Tadmor 
18fl is based on the third-order reconstruction by Liu and Osher in | i~6|| . For the two- 



dimensional method see fjj and 0. 

A first step for importing the high-order reconstructions that were derived in the 
upwind framework was taken in 0. There, the ENO method was transformed into 
the central setup, and a new mostly centered stencil was shown to produce the least 
oscillatory results. The next step was taken in the ID case in |T2|], where a new cen- 
tral weighted non-oscillatory (CWENO) reconstruction was introduced. This CWENO 
method is based on the upwind WENO methods by |17| and ||, in which an inter- 
polant is written as a convex combination of several reconstructions which are based 
on different stencils. These methods include an internal switch which is designed such 
as to provide the maximum possible accuracy in smooth regions, while automatically 
switching to a the more robust one-sided stencil in the presence of discontinuities and 



large-gradients. The ID CWENO method was extended to two space dimensions in |L4 



and ||15||. For a numerical study of the behavior of the total variation for the CWENO 
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scheme, see fO 



In this paper we present a new, compact CWENO reconstruction. This new re- 
construction is based on defining an arbitrary quadratic function which is added to 
linear interpolants in such a way as to obtain third-order accuracy in smooth regions 
(in one and two space dimensions). In regions with discontinuities or large gradients, 
the weights are automatically changed so that they switch to a one-sided second-order 
linear reconstruction. These reconstructions turn to be extremely compact; in the one 
dimensional case, e.g., the reconstruction is based on a four-point stencil. 

The structure of this paper is as follows: We start in §|2| with a brief overview of 
central schemes for conservation laws in one and two space dimensions. 

We then proceed to present our new, compact, third-order CWENO reconstruction 
in §£| The idea of introducing an arbitrary quadratic reconstruction is first presented 
in the one dimensional framework and then extended to two space dimensions. 

We end in §|3] where we demonstrate our new method in several test cases. First, 
the accuracy tests (both in one-dimensional and two-dimensional cases) show the third- 
order accuracy of the method. We then solve the one-dimensional system of the Euler 
equations of gas dynamics for a few test problems and we illustrate the behavior of 
the scheme in scalar two-dimensional cases. In particular, we would like to stress that 
in our numerical results we observe a very robust and non-oscillatory behavior of the 
weights, which can be related to the overall robustness and accuracy properties of our 
new method. 
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under contract DE-AC03-76-SF00098. Part of this work was done while CP. and G.R. 
were visiting the Lawrence Berkeley Lab. 



2 Central Schemes for Conservation Laws 

In this section we give a brief overview of central schemes for approximating solutions 
to hyperbolic conservation laws in one and two space dimensions. For further details we 
refer the reader to [23], M, [12| and the references therein. 

Starting in the one-dimensional case, we seek numerical solutions of the Cauchy 
problem 

u t + f(u) x = 0, 

(2.1) 

u(x, t = 0) = uo(x). 

For simplicity, we introduce a uniformly spaced grid in the (x, t) space, where the mesh 
spacings are denoted by h := Ax and k := At, respectively. We denote by w™, the 
numerical approximation of the cell average in the cell Ij := [xj-i/2, Xj+1/2] at time 
t n = nk, where Xj = jh. The finite-difference method will approximate the cell-averages 
at time t n+l based on their values at time t n . 
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We start by reconstructing at time t n a piecewise-polynomial conservative interpolant 
from the known cell- averages, uj , i.e., 

3 

where Xj is the characteristic function of the interval Ij and Rj(x) is a polynomial de- 
fined in Ij. It is here, in the reconstruction step, where the accuracy and non-oscillatory 
requirements enter. Different central methods will be typically based on different recon- 
structions. 

The reconstruction, P u (x, t n ), is then evolved exactly in time (integrating (|2.ip), and 



then projected on staggered cells, in order to compute the cell average at Ij+i/2- With 
this procedure we obtain 



i7 n+1 - il n 
U j+l/2 - "j+1/2 



i r +1 

+ - J \f{Pu{xj, r)) - f(P u (x J+1 , t))} dr. (2.3) 



Based on the reconstruction, Q2.2|) , the first term on the RHS of (|2.3|) can be explicitly 
computed, 



I r x 3+i 

^j+l/2 = T I P u (x,t n )dx. 

J Xj 



In order to compute the second integral on the RHS of ( |2.3p , namely the integral in time 
over the fluxes, one should observe that due to the staggering, up to a suitable CFL 
condition, these integrals involve only smooth functions, and hence can be approximated 
by a sufficiently smooth quadrature. A second-order method can be obtained, e.g., using 
the mid-point rule in time (see ||19|| ). Simpson's rule for the quadrature in time will 
provide fourth order accuracy (which will naturally be sufficient also for a third-order 
method). 

The quadratures for the time integrals of the fluxes, require the prediction of point- 
values of the function in several points in the interval [t n , t n+1 ]. One possible approach is 
to use a Taylor expansion based on the equation, ( |2.1| ). Such an approach was used, e.g., 
in and fl8f . In order to avoid the technical complications involved in the Taylor series 



(in particular with high-order methods, and when dealing with systems of equations), 
one can alternatively use a Runge-Kutta (RK) method directly on ( |2.1| ), to predict the 
required values in later times. By using the Natural Continuous Extension (NCE) of 



Runge-Kutta schemes [£5| , the value of the flux at the nodes of the quadrature formula 
can be computed with a single RK step. Such a method is simpler compared with the 
Taylor expansion method, but it does require another reconstruction: the reconstruction 
of the point values of the derivatives of the fluxes at time t n which are then used as the 
input to the RK solver (see and |T2j for more details). 



The simplicity of central schemes manifests itself when turning to deal with systems 
of equations. Basically, the algorithm that was described in the scalar case, repeats 
itself component-wise. Based on the type of the reconstruction, when solving systems 
of equations, there can even be simplifications over a purely componentwise extension 
of the scalar scheme. These can be found, e.g., in §|3~3| below. 
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The extension to two space dimensions is straightforward; We consider 

u t + f(u) x + g(u) y = 0, (2.4) 

subject to the initial condition, u(x, y, t — 0) = Uo(x,y). Here, Ax and Ay will denote 
the spatial mesh spacings, while At denotes the spacing in time. The two-dimensional 
cells are now I ifj = [0^-1/2,^+1/2] x [%-_i/ 2 , Jfr+i/a]- 

Following the general methodology, we wish to construct the cell-averages, u^ 1 , 
based on the cell averages at time t n . First, we construct a two-dimensional interpolant 
which reads 

P u (x, y, t n ) = v)Xi,j, (2-5) 

with Xi,j being the characteristic function of the cell Iij and Rij(x,y) a polynomial of 
a suitable degree. An exact evolution in time of the interpolant fl2.5|) which is projected 
on its cell-averages now reads (compare with (|2.3|)) 



U i+l/2,j+l/2 ~ M r+l/2,j+l/2 + (2-6) 

+ a„a- / / [f(Pu(xi,y,T)) - f(P u (x i+1 ,y,r)]dydr + 
Jy=yj 



Ax Ay J t n 



+ A a / / [f(Pu{x, yj, r)) - f(P u {x, y j+ i, t)) dxdr. 

LSXLiy J t n J x =Xi 

The staggered cell-average at time t n can be directly computed by 



AxAy 



< + i/ 2 ,j+i/2 = a_ a .. / / P u (x,y,t n )dydx. 



Analogously to the one-dimensional case, the flux integrals on the RHS of (pTq) can be 
approximated using a quadrature rule in time. This can be coupled with a quadrature 
rule for the line integrals in space. Here it is necessary to avoid quadrature points at 
which the fluxes may not be smooth. The details can be found, e.g., in [|14], [15| (see also 



3 A Compact third-order CWENO Reconstruction 
3.1 The One-Dimensional Framework 

In this section we derive our new CWENO reconstruction in one space dimension. The 



two dimensional extension will follow in §3.2 



We first note that in the absence of large gradients, we obtain third order accuracy 
if we choose for the reconstruction the optimal polynomial 
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where P OPT j(x) is the parabola that interpolates the data tt"_ 1; u", in the sense of 
cell averages to enforce conservation: 



x j+l+l/2 

Poptj(x) dx = u] +l , I = -1, 0, 1. 
These conditions determine P OVT j(x) completely, namely: 

Po*tA x ) = u ] + u 'j( x - x j) + \ u "M - x if> 

with 

«" = fi?-^(^i-2^ + fl?-i), 

J ' 2Ax ' J ' Ax 2 

However, when discontinuities or large gradients occur, this reconstruction would be 
oscillatory. Therefore following the WENO methodology (|17j) @, ||12||)) we construct 



an essentially non-oscillatory interpolant as a convex combination of polynomials which 
are based on different stencils. Specifically, in the cell Ij we write 

!)(.,■) ^^/y(.r). $>f = l, Wi >0, z6{ l ,c,k}, (3.2) 

i i 

where P L and P R are linear functions based on a left stencil and a right stencil, respec- 
tively, and P c is a quadratic polynomial. In order to simplify the notations, we will omit 
the upper index j, remembering that the weights and the three polynomials change from 
cell to cell. 

Conservation requires that P R (x) will interpolate the cell averages 



/ P R (x)dx — AxuJ, / P R (x)dx = AxUj +1 , 



'Xj-1/2 Jx ] + l/2 

which in turn, results with 



Pn(x) = u] + -^r J -(x - Xj ). (3.3) 

J Ax 

Similarly, for the left interpolant we have 



P L (x) = u] + 3 Ax J ( x - x j)- (3-4) 

All that is left is to reconstruct a centered polynomial, P c , such that the convex combi- 
nation, (|3.2|) , will be third-order accurate in smooth regions. It must, therefore, satisfy 



P OVT (x) = C h P^x)+C R P R {x)+C c P c {x), ^Ci=l, z€{l,c,r}, (3.5) 

i 

where C h , C c and C R are constants. Due to the staggering between every two consecutive 
steps of the central method, our reconstruction should provide half-cell averages which 
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are third-order accurate. A straightforward calculation shows that any symmetric choice 
of constants in ( |3.5|) provides the desired accuracy. In particular, for the specific 
choice of C L = C n = 1/4, equations (^j.3|) (PT5|) yield 

P c (x) = 2P OPT (x)-^(P R (x) + P L (x))=u]-^(u] +1 -2u] + u]_ 1 ) + 

+ ' 2Ax ~ ^ + AX 2 — (* - • ( 3 - 6 ) 

In order to complete the reconstruction of Pj(x) in (|3~^ ), it is left to compute the 
weights Wi. Following |J, |12[, we write 

Wi = v , Q!l , Q< = 7 — j= v- , e {l,c,r}. (3.7) 



The constants CVs in ( j3/7|) are chosen to be the same as in (|3.5|) , i.e., C L = C R = 1/4, 
while C c = 1/2. 

The smoothness indicators , ISj, are responsible for detecting large gradients or dis- 
continuities and to automatically switch to the stencil that generates the least oscillatory 



reconstruction in such cases. Once again, we follow |J, [|1J] and define in each cell Ij 
the three smoothness indicators, ISj, as 



+V2 h^iP^ix^dx. % E {l, c, a}. (3.i 



A direct computation, based on (p.3|), (|3.4|) and (|3.6|) , yields 

IS L = (uj-u]^) 2 , IS R = (u] +1 -u])\ 

1 3 1 

ISa = Y (u] +1 -2u] + u]_ 1 ) 2 + -(u] +1 -u]_ 1 ) 2 . (3.9) 

The constant e is taken as to prevent the denominator from vanishing. Furthermore, its 
valued has to satisfy two requirements, namely 

i) e ^> IS in smooth regions 

ii) e <C IS near discontinuities 

The first conditions guarantees that, on smooth regions, the weights are basically equal 
to the constants that provide high accuracy. The second conditions guarantees that in 
the presence of a discontinuity, the weights of the parabola and of one of the one-sided 
linear reconstructions will be practically zero, and we will be left with the other one- 
sided linear reconstruction. Hence, our third-order method automatically switches to a 
second-order method in the presence of large gradients, which is exactly what makes it 
so robust as will be evident in our numerical computations presented below. 

The constant p weights the departure from smoothness. We used the value p = 2. 
For a discussion about its choice see || and |T^| . 

A second, non-oscillatory reconstruction is required for the flux derivative. It is 
natural to adapt the reconstruction that we used for the half-cell averages also to the 
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reconstruction of the point-values of the flux derivative. Here however the interpolation 
requirements will be in the sense of point values instead of cell averages. Once again, 
a direct computation shows that any symmetric choice of constants will provide the 
desired accuracy and it will only be natural to use the same constants, Cj's, that were 
used in ( |3.5| ). This is simpler than the case of [|EJ where we had to use different sets 
of constants for the two different reconstructions (the reconstruction of the half cell 
averages, and the reconstruction of the point-values of the flux derivative at the edges 
of the domain). 



Remarks: 



We would like to emphasize that our new method is based on adding the arbitrary 
parabola P c into the convex combination which is the heart of our non-oscillatory 
reconstruction. Our numerical simulations showed that the freedom we have in 
selecting the constants Cj has no influence on the properties of the method. It is 
easy to prove that we obtain third-order accuracy regardless of the smoothness of 
the weights, as long as they are symmetric. This is substantially more robust than 
the third-order method of Liu and Tadmor in where the order of the method 



did depend on the smoothness of the limiters, and could deteriorate to first order 
in the presence of large gradients (as was shown in @). 

By extending our new ideas, one can modify our previous CWENO method pre- 
sented in [HJ in order to obtain a fifth-order, central, non-oscillatory scheme. This 



can be done by simply adding a fourth-order polynomial for the computation of 
the point-values. 

3. Our new reconstruction is equivalent to limiting with the minimum slope instead 
of slope zero in the presence of a discontinuity. Hence, it is based on a second-order 
reconstruction close to shocks, unlike the scheme of [ITBJ which can be only first 
order accurate in such regions. 

4. In the one dimensional case, the additional parabola is needed only for the accu- 
rate recovery of the point-values. In the 2D framework, however, the equivalent 
additional parabola will be required also for the accurate reconstruction of the 
fractional cell-averages. 



3.2 A Two-Dimensional Extension 



We extend the ideas of § |3.1| to the two-dimensional framework. This extension is 
straightforward and is based on reconstructing an interpolant as a convex combination 
of four one-sided, piecewise- linear interpolants, and a centered, quadratic interpolant, 
such as to get the desired third-order accuracy in smooth regions. 

Following these ideas, the reconstruction in the cell Jy, can be written as 



Pi,j(x,y) = ^W^P^ix^), k E {NE,NW,SE,SW,c}, 

k 



(3.10) 
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with ^2 k w^ J = 1, and w^ J > 0. Here, P NE , P NW , P SE and P sw are the one-sided linear 
reconstructions, and P c is a centered quadratic reconstruction (see Figure |3.1|) . Similar 
to the one- dimensional case, we will simplify our notations by omitting the superscripts, 
remembering that both the weights and the polynomials, change from cell to cell. 



NW 








NE 








/ 




















sw 








SE 



Figure 3.1: The Two-Dimensional Stencil 

Clearly, in an analog to the one-dimensional case, ( |3.3| - |3.4|) , the four linear recon- 
structions are given by 

P NE (x,y) = u id + ^ — -(x-Xi) + A ^ — -(y-yj), 

P«w(x,y) = <, • ' A ; Ll (■'■■-,,) • -' (y .///). 

P*«(x,v) = ul 3 + ^ ~ ' ' (x - x.) + ^ (y - yi ), 

= <■+ ^ ^-^± {y- yj ). (3.11) 

The centered polynomial, P c (x,y), is taken such as to satisfy 

-Popt(^, 2/) = C kPk(x, y), ) j C k = l, k e {ne, nw, SE, sw, c}. (3.12) 
k k 

Here, P pt is the quadratic polynomial based on a nine-point stencil, centered around 
Jjj, which is given by (see |p"T| ) 

Popt (x,y) = Ui d + u' itj (x-Xi)+ u\ j (y - y, ) + u'-j ( x - x t ) (y - yj ) + 

+^ul j (x-x i ) 2 + ±ul j (y-y j ) 2 , (3.13) 

where 



= - i ((Ax) a i^ + (Ay) a «y , 

■7/"' 1 1^" 11^ II 
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u 



Ax 2 

n I K-n 

i+lj+l + U i-X,j-1 



Ay 2 



(3.14) 



U i+l,j-l U i-l,j+l 



AAxAy 

Unlike the one-dimensional case, not every symmetric selection of the constants CVs 
will provide a third-order reconstruction for the quarter cell- averages. Here, a straight- 
forward computation shows that in order to satisfy the accuracy requirements, we must 
take C NE = C NW = C sw = C SE = 1/8. Hence, C c = 1/2, and (|3.12|) implies 



P c (x,y) = 2P Q 



:(x, y) - - P NE {x, y) + P NW (x, y) + P sw {x, y) + P SE (x, y) 



+o u "A x ~ x i) 2 + o u lj(y - %) 2 ' 



(3.15) 



where 







- (A 
12 L v 






u h = 


II 


= 2 <i> 


u h = 



{Axful 3 + {Ayfu 



2„.w 



- 2 <p 



u 



h3 



2 <y 



All that remains is to determine the weights w k ,J in ( |3.10 ). Once again, we write 



w 



a 



Li a i v- -T- ^ k 

The constants, Ck, are the same constants that were used to reconstruct the centered 
parabola in ( ft. 12 ). The constants, e and p, play the same role as in the one- dimensional 
case. At that point, to simplify the notations we assume that the mesh spacings are 
equal in the x and y directions, i.e., Ax = Ay = h. We can then follow [14], and define 
the smoothness indicators, IS]f , as 

IS ^ = / / h 2 ^~ 1 \D a P k ) 2 , ke {ne,nw,se,sw,c}. (3.16) 

|a|=l,2 ^ x i-h/2 Jyj-h/2 

If Aa; Ay, then only a trivial enhancement to ( ft,16| ) is required. For the four one-sided 
linear reconstructions, which can be all written as 

Pk = u + u(x — Xi) + u(y — yj), k e {ne, nw, se, sw}. 

with suitable reconstructed point-values and first derivatives, a direct computation of 
( ft. 161 ) results with 

IS k = h\u'f + {u v ) 2 }. (3.17) 

The centered smoothness indicator, IS C , which corresponds to the centered quadratic 
reconstruction, P c (x,y), ( ft. 151 ), is given by 

IS C = h 2 [{u'f + {u x ) 2 ] + %- [13( M ") 2 + 14K) 2 + 13K) 2 ] , 

o 

(the discrete derivatives are given by ( ft,15| )). 



(e + ISl' 3 > 



k,le{ 



NE, NW, SE, SW, C 
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3.3 A Note on Systems 

Almost nothing changes when turning to deal with systems. The reconstruction that 
was described in the previous sections, directly transforms to systems and is performed 
component-wise. The only relatively subtle issue is the computation of the smoothness 
indicators. 



In our previous work we have suggested several approaches for the computation 
of the smoothness indicators. Three different options were suggested: the first is to allow 
every component to have a strictly individual behavior, namely to allow a different stencil 
with different smoothness indicators to each component. In the second approach, we 
designed global smoothness indicators such as to force every component to adjust even 
when the discontinuity is in a different component. The last approach was to use external 
information about the system. For example, in the Euler equations of gas dynamics, 
one expects both shocks and contact discontinuities to occur in the density. Hence, all 
stencils can be tuned according to this component. 



Our results in [|1^] showed that the best approach is the one which was based on 
the global smoothness indicators. It requires no additional information on the system, 
it produced less oscillatory results compared with the individual smoothness indicators 
for each component, and it was the simplest to implement. 

In the one-dimensional case, e.g., the global smoothness indicators are given by 
(compare with 



is k = i Yl fit ( E f 3+1/2 h2l ' x ( p £) 2 dx i ' k G k c ' R > ( 3 - 18 ) 

Here the k-th polynomial in the r-th component is denoted by Pk, r , and d is the number 
of equations. The scaling factor \\u r \\2 is defined as the L 2 norm of the cell averages of 
the r-th component of u, 




The numerical examples performed for systems, appearing in §|4], were carried out 
with the global smoothness indicators. 



4 Examples 



We present numerical tests in one and two space dimensions, in which we demonstrate 
the accuracy, robustness and high-resolution properties of our new method. 

|14|| and in all our numerical examples 



Following our previous works 



12 



we integrate in time using a Runge-Kutta method with natural continuous extension, 



which was presented in p5[, with a fixed time step. 
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The time step is determined by imposing that the Courant number is a given fraction 
of the maximum Courant number determined by linear stability analysis. The Courant 
number is defined by 

c = ^ 

maxj pj 

where pj denotes the spectral radius of the matrix f'(v,j) computed on the initial con- 
dition, and A = At/ Ax is the mesh ratio. 

The linear stability analysis carried out in [[| yields a Courant number C = 3/7 
for the one dimensional case. We remark that the stencil used in M was different than 
the one that we use, and therefore linear stability should be repeated for the compact 
scheme in order to obtain a sharp estimate on the maximum Courant number. 



Example 1: Accuracy Tests 

Our first example checks the accuracy of our new method in several one and two- 
dimensional test cases. In all of the one-dimensional tables, the norms of the errors are 
given by 



L 1 — error : 1 1 Error 1 1 1 = YljLi \ u ( x j>t n ) ~ u^\h 



L°° — error : 1 1 Error = maxi<j<jv \u(Xj,t n ) — u™\. 
Analogous expressions hold for the two-dimensional norms. 

1. Linear advection: This test estimates the convergence rate at large times. We 
solve u t + u x = 0, subject to the initial data u(x,t = 0) = sin(7rx) and to periodic 
boundary conditions on [—1, 1]. The integration time was taken as T = 10. 

In Table |4.1| , we show the results obtained for this test problem with e = 10~ 2 . 
We clearly see a third order convergence rate in both L l and L°° norms. We also 
show in Table |4.2| the results obtained for the same example with e = 10~ 6 , which 



is the value suggested in both i| and [JI2" 

Compared with Table the errors are larger, and the convergence rate is not as 
regular as before. This is mainly due to the fact that for a very small value of s, 
condition i) is not satisfied, until the grid spacing Ax becomes very small. 

2. Linear advection with oscillatory data: This test is used to detect deteriora- 
tions of accuracy due to oscillations in the parameters that control the selection of 
the stencil (for details see || and the references therein). Once again, the equation 
is u t + u x = 0, subject to the oscillatory initial data, u(x, t — 0) = sin 4 (7ra;) and 
to periodic boundary conditions on [—1, 1]. Here, the integration time is taken as 
T = 1 and e = 1(T 2 . 



The results of this test are displayed in Table 4~3[, and confirm the third-order 
accuracy of the method with no deteriorations in its accuracy. 
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3. Burgers equation: We solve the Burgers equation u t + (0.5m 2 ) z = 0, subject to 
the initial data u(x,t = 0) = 1 + 0.5sin(7rx) and to periodic boundary conditions 
on [—1, 1]. The integration time is T = .33, and e = 1CT 2 . Here, a shock develops 
at T = 2/tt. Note that here the maximum speed of propagation is f'(u) = 3/2. 
Thus we use A = .66 * 3/7 ~ 2/3A max 



Table [4.4| shows the results we obtained which verify the third-order accuracy of 



the method also for nonlinear problems. 

4. 2D Linear advection: Finally, we implemented our method for the two-dimensional 
linear advection problem, u t + u x + u y = 0. The initial condition is taken as 
u(x, t = 0) = sin 2 (7rx) sin 2 (7n/), and we impose periodic boundary conditions on 
[0, l] 2 . The errors and the estimated convergence rate our computed at time T = 1. 



In Table fig we present the results obtained when the weights are taken as con- 
stants ( |3.12| ). Table [4. 6| shows the results obtained with the fully non-linear 



scheme, where the weights of the reconstruction include also the oscillatory in- 
dicators. With constant weights our method is third-order as expected, while 
with non-constant weights, there seems to be a better convergence rate. A careful 
study of the tables shows, however, that this better convergence rate is mainly 
due to larger errors on the coarser grids. 



N 


ji error 


L 1 order 


L°° error 


L°° order 


20 


0.1423 




0.1484 




40 


0.1308E-01 


3.44 


0.1708E-01 


3.12 


80 


0.7054E-03 


4.21 


0.1071E-02 


4.00 


160 


0.7517E-04 


3.23 


0.7823E-04 


3.78 


320 


0.9391E-05 


3.00 


0.7977E-05 


3.29 


640 


0.1174E-05 


3.00 


0.9406E-06 


3.08 


1280 


0.1467E-06 


3.00 


0.1158E-06 


3.02 



Table 4.1: Linear advection, T = 10, Uq(x) = sin(7nc), e = 10 2 , Courant number 
C = 0.9C max , C max = 3/7 



Example 2: ID Systems - Euler Equations of Gas Dynamics 

Next, we solve the Euler equations of gas dynamics, 




rn 



pu 2 + p =0, p=( T -l). (e-^-u 2 ) 
u(E + p) 1 - 



where the variables p, u, m = pu, p and E denote the density, velocity, momentum, 
pressure and total energy, respectively. We use two sets of initial data: 
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1\T 
IN 


error 


L 1 order 


L°° error 


L°° order 


20 


0.2292 




0.2522 




40 


0.8975E-01 


1.35 


0.9943E-01 


1.34 


80 


0.2184E-01 


2.04 


0.3759E-01 


1.40 


160 


0.3677E-02 


2.57 


0.1090E-01 


1.79 


320 


0.3682E-03 


3.32 


0.1896E-02 


2.52 


640 


0.2454E-04 


3.91 


0.1585E-03 


3.58 


1280 


0.1379E-05 


4.15 


0.5972E-05 


4.73 



Table 4.2: Linear advection, T = 10, Uq(x) = sin(7rx), e = 10 6 , Courant number 
C = 0.9C max , C max = 3/7 



N 


j\ error 


L 1 order 


L°° error 


L°° order 


20 


0.1285 




0.1909 




40 


0.2813E-01 


2.19 


0.5223E-01 


1.87 


80 


0.2608E-02 


3.43 


0.5226E-02 


3.32 


160 


0.2553E-03 


3.35 


0.3619E-03 


3.85 


320 


0.3055E-04 


3.06 


0.3319E-04 


3.45 


640 


0.3826E-05 


3.00 


0.3814E-05 


3.12 


1280 


0.4777E-06 


3.00 


0.4654E-06 


3.03 



Table 4.3: Linear advection, T — 1, Uo(x) = sin 4 (7r:r), e = 10 2 , Courant number 
C = 0.9C max , C max = 3/7 



N 


L 1 error 


L 1 order 


L°° error 


L°° order 


20 


0.7974E-02 




0.1527E-01 




40 


0.6654E-03 


3.58 


0.1844E-02 


3.05 


80 


0.6563E-04 


3.34 


0.2340E-03 


2.98 


160 


0.8494E-05 


2.95 


0.3645E-04 


2.68 


320 


0.1067E-05 


2.99 


0.4937E-05 


2.88 


640 


0.1355E-06 


2.98 


0.6388E-06 


2.95 


1280 


0.1695E-07 


3.00 


0.8047E-07 


2.99 J 



Table 4.4: Burgers equation, T = .33, u (x) = 1 + 0.5 sin(7ra;), e = 10~ 2 , A = 0.66 * 3/7, 
corresponding to a Counant number C = 0.99C max 
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N 


ji error 


L 1 order 


L°° error 


L°° order 


10 


3.696E-02 




1.252E-01 




20 


4.964E-03 


2.90 


1.767E-02 


2.83 


40 


6.304E-04 


2.98 


2.264E-03 


2.96 


80 


7.902E-05 


3.00 


2.842E-04 


2.99 


160 


9.880E-06 


3.00 


3.555E-05 


3.00 



Table 4.5: 2D Linear Advection, Constant weights; T = 1, A = 0.425, uq(x) = 
sin 2 (7ra) sm 2 (ny), 



N 


L l error 


L 1 order 


L°° error 


L°° order 


10 


9.750E-02 




3.447E-01 




20 


1.419E-02 


2.78 


8.111E-02 


2.09 


40 


9.387E-04 


3.92 


7.967E-03 


3.35 


80 


8.319E-05 


3.50 


4.465E-04 


4.16 


160 


9.977E-06 


3.06 


3.999E-05 


3.48 



Table 4.6: 2D Linear Advection, Non-Constant weights; T = 1, A = 0.425, Uq(x) = 
sin 2 (7rx) sm 2 (iiy), e = 10~ 2 



1. Shock tube problem with Sod's initial data, [23 



{p l ,m h E l ) = (1,0,2.5), x<0.5, 
(p r , m r , E r ) = (0.125, 0, 0.25), x > 0.5. 



2. Shock tube problem with Lax initial data, 



{Pr 



(0.445,0.311,8.928), 
= (0.5,0,1.4275), 



x < 0.5, 
x > 0.5. 



We integrate the equations up to time T = 0.16. In Figure |4.1| we show the density 
components for Sod's initial data, and in Figure |4.2j we show the equivalent plot for 
Lax initial data. In both figures we also present the weight of the central stencil at the 
final time. All the results are given for 200 and 400 grid points. Following |19| , we pick 
A = .1. Note that the maximum characteristic speed for Sod's problem is roughly 2.5, 
while for Lax problem the maximum propagation speed is ~ 5. Since our scheme has a 
Courant number no larger than .5, we see that A = .1 is actually the maximum value 
compatible with stability. This explains while there are still some wiggles in the test 
solution of Lax, while the Sod's solution is monotone. 

It is interesting to compare the behavior of the central weight of the new method to 
the behavior of the central weight of the original CWENO method \T2\. Here, the weights 
are much smoother compared with the weights in [12]. The accuracy and stability 



properties of the method can be related to the smoothness of the nonlinear weights 
involved (see [0]). 
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Density, N=200 



Density, N=400 




0.2 0.4 0.6 0.8 1 




0.2 0.4 0.6 0.8 1 



0.8 
0.6 
0.4 
0.2 



Central Weight, N=200 




0.8 
0.6 
0.4 
0.2 



Central Weight, N=400 



V 



0.2 0.4 0.6 0.8 1 



0.2 0.4 0.6 0.8 1 



Figure 4.1: Euler equations of gas dynamics - Sod initial data, A = 0.1, T = 0.16 



Density, N=200 



Density, N=400 





0.2 0.4 0.6 0.8 1 



0.2 0.4 0.6 0.8 1 



Central Weight, N=200 



Central Weight, N=400 





0.2 0.4 0.6 0.8 1 



Figure 4.2: Euler equations of gas dynamics - Lax initial data, A = 0.1, T = 0.16 
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Example 3: 2D Problems 

1. Linear rotation: Following fC4]| , we consider a linear rotation of a square patch on 
[0, l] 2 , with initial condition 1*0(2;, y) = 1 for {\x- 1/2| < 1/2} x {|y- 1/2 1 < 1/2} 
and zero elsewhere. In Figure [Q| we display the solution after a rotation of 7r/4 and 
of 7r/2. There are no spurious oscillations. We also show the corresponding central 
weight. As expected, this weight is zero in regions where the solution has steep 
gradients, and that is exactly the property that prevents spurious oscillations from 
developing. Even though we are dealing with linear waves, the resulting resolution 
is relatively good. Due to the compactness of the stencil, when the slopes are not 
sharp, they are not identified as discontinuities. This can be observed in the plots 
of the central weight, which returns to its constant value (1/2) on the slope. 





Figure 4.3: Linear Rotation, A = 0.425, iV = 40 



2. 2D - Burgers equation: We end by solving the two-dimensional Burgers equa- 
tion, u t +(u 2 /2) x +(u 2 /2) y = 0, subject to the initial dataw(x,t = 0) = sin 2 (7rx) sin 2 (7n/) 
and periodic boundary conditions on [0, l] 2 . In Figure [Ojwe present the solution 
obtained at time T = 1.5 with different mesh spacings. One can easily notice that 
the shocks are well resolved and there are no spurious oscillations. 
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Figure 4.4: Two-Dimensional Burgers equation, A = 0.425, T = 1.5 
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